Vortex nucleation in mesoscopic Bose superfluid and breaking of the parity symmetry 
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We analyze vortex nucleation in mezoscopic 2D Bose superfluid in a rotating trap. We explicitly include a 
weakly anisotropic stirring potential, breaking thus explicitly the axial symmetry. As the rotation frequency 
passes the critical value £l c the system undergoes an extra symmetry change/breaking. Well below Q c the 
ground state is properly described by the mean field theory with an even condensate wave function. Well above 
Q c the MF solution works also well, but the order parameter becomes odd. This phenomenon involves therefore 
a discrete parity symmetry breaking. In the critical region the MF solutions exhibit dynamical instability. The 
true many body state is a strongly correlated entangled state involving two macroscopically occupied modes 
(eigenstates of the single particle density operator). We characterize this state in various aspects: i) the eligibility 
for adiabatic evolution; ii) its analytical approximation given by the maximally entangled combination of two 
single modes; and finally iii) its appearance in particle detection measurements. 



PACS numbers: 03.75.Hh, 03.75.Kk, 67.40.Vs 

I. INTRODUCTION 

Symmetry changes or breaking belong to some of the most 
fascinating phenomena in nature. In classical physics they are 
often associated with phase transitions in macroscopic sys- 
tems (DEI. Paradigm examples of symmetry breaking con- 
cern magnetic phenomena, such as for instance appearance of 
ferromagnets at temeperatures lower than the, Curie tempera- 
ture, T c . In the classical world symmetry changes/breaking 
(C/B) are driven by thermal fluctuations, and in the stan- 
dard Landau-Ginsburg scenario are associated with increase 
of classical correlations and arousal of the long range order. 
Mean field approach, that goes back to "molecular field the- 
ory" of Curie-Weiss |3|. provides very often quite correct de- 
scription of these phenomena away from criticality. Close to 
critical temperature, quantitative description requires the use 
of renormalization group approach a la Wilson |H1[5). 

In quantum physics paradigm examples of symmetry C/B 
deal with low temperature behavior of weakly interacting 
quantum Bose gases and Bose-Einstein condensation (BEC) 
(6). In the quantum world particularly interesting are quan- 
tum phase transitions [7|, and quantum symmetry C/B that 
are driven by quantum fluctuations. They can occur either at 
zero temperature, or in quantum dynamical externally driven 
systems. 

The symmetry C/B that have drawn a lot of attention since 
the early discovery of superfluids |8| till the recent studies 
of BEC is nucleation of vortices in rotating superfluids. In 
fact, one of the most striking properties of superfluid and con- 
densed systems is their response to rotation. The only way to 
acquire angular momentum is by the nucleation of vortices, 
topological singularities surrounded by condensed atoms re- 
volving around their cores. The cores are well localized and 



have size of the, so called, healing length. At low temperatures 
they are empty, whereas at higher temperatures they are filled 
with the thermal fraction of the condensate. For quantum 
gases, atoms are usually confined in an isotropic harmonic 
trap and experience an additional quadratic potential rotating 
at angular frequency (for a review see ||9] [TO)). Standard 
textbooks |6 1 associate vortex nucleation with thermodynamic 
instability. When the rotation frequency is small, there exists 
a mean field solution of the equation describing the BEC or- 
der parameter (condensate wave function) with a single vortex 
lfTTlfT2 l: this solution has, however, larger (free) energy than 
the one corresponding to the condensate at rest fl3l . Above 
certain critical rotation frequency, the solution with the vortex 
becomes a ground state, and in principle may be achieved at 
low temperature being driven by thermal fluctuations. In prac- 
tice, experiments with BEC occur in a completely different 
way. Typically, one prepares a condensate at very low tem- 
perature, and then applies a certain dynamical perturbation to 
create a vortex. First vortices have been created at JILA lfl4l 
using a kind of phase imprinting method ifTSIl . It turned out, 
however, that the method consisting of slight deformation and 
rotation of the trap, or alternatively, "laser stirring" |[T6l was 
more efficient and led to numerous spectacular observations 
such as that of Abrikosov lattice (|[T6l[T7l[T8l[T9ll20ll2T1l22l. 
for a review see l9l l23l ). 

Bose-Einstein condensates (BEC's) of dilute atomic gases 
offer particular possibilities for studies of nucleation of 
vortex-states and their expansion in the course of time-of- 
flight (TOF) detection. In addition, these system, allow for 
experimental analysis and manipulation of mesoscopic con- 
fined clouds of condensed atoms trapped in the sites of op- 
tical lattices. This opens the perspective to compare directly 
results of exact numerical analysis for small systems with ex- 
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periments and with different approximate calculations (for the 
first experiments in this direction see 11241 '). Particulalry inter- 
esting are the studies of the applicability of the mean field 
(MF) theory within specific conditions. For confined systems 
at large rotation frequencies (in analogy to charged particles 
submitted to high magnetic fields) strongly correlated states 
develop, the so called fractional quantum Hall (FQH) type 
states flOl l25l l26l . Such states cannot be described by a sin- 
gle particle wave function that would play the role of an order 
parameter. 

From a theoretical point of view, the vortex nucleation can 
be tackled by several techniques, ranging from a MF approach 
based on the Gross-Pitaevskii (GP) equation E71 l28l 1291 l30l 
to the investigation of the many-body energy eigenstates 
EU|32l|3l|34ll35l|36l|33. Many authors addressed the ques- 
tion of vortex nucleation theoretically, asking in the first place 
for energetic stability of the vortex configuration as a ground 
state. Within the mean field approach this has been discussed 
mainly in the context of thermodynamic stability (cf. ll6l [T3"ll ). 
Several papers discussed, however, the case of T = and 
vortex nucleation in the ground state (GS) of the system using 
the exact quantum description (cf. (25] [32] [34j [35] [38]), or 
rigorous derivation of the MF equations [39, 40 1. More recent 
papers treated the problem of dynamics of vortex, or vortex 
lattice nucleation in elliptically deformed rotating traps, us- 
ing mean field method (i.e., time dependent Gross-Pitaevski 
equation (GPE) [6]) and trying to reproduce the experimental 
results. The conclusion of these works is that vortex nucle- 
ation is inevitably associated to dynamical instability of the 
solutions of GPE l28l 141) . Same results hold in the case of 
vortex nucleation via phase imprinting |42|. Some authors 
||29l claim that apart from the dynamical instability, a Landau 
instability (associated to dissipation) is also necessary to allow 
vortices to penetrate the BEC. The dynamical instability GPE 
is genetically associated to the appearance of squeezing of two 
Bogoliubov-de Gennes (BdG) quasi-particle modes, i.e., ex- 
ponentially growing two mode entanglement in the regime of 
validity of BdG (i.e., regime of small Gaussian fluctuations 
around the MF solutions) P3"]| . This observation already in- 
dicates the necessity of going beyond the mean field at the 
instability. 

Much less is known about exact dynamics of the vortex 
nucleation. Parke et al. [37 1 considered recently this prob- 
lem for a mesoscopic sample of atoms in the lowest Landau 
level (LLL) and discover striking non-mean field effects in 
the (stationary) spectrum of the system at the critical rotation 
frequency Q c . They interpreted their results in terms of coop- 
erative tunnelling of a vortex pair by "requantizing" the mean 
field theory (reduced to 3 relevant modes). In a recent paper 
PRl we used another approach and studied exact dynamics of 
a mesocopic sample of atoms in a elliptically deformed rotat- 
ing harmonic trap. Our main result was that as one increases 
f2, at the f2 c the mean field description ceases to be valid. The 
system enters a strongly correlated and entangled state, well 
described by an effective two mode model. The mean field 
description (similar to that of ref. Il3~7"ll ) exhibits dynamical 



instability and hysteresis for il ~ fl c . Since we explicitly in- 
clude an anisotropic stirring potential, the present mechanism 
concerns a discrete parity symmetry breaking. Therefore it 
differs from the case of the vortex nucleation in axially sym- 
metric traps: in the latter case, breaking of the continuous ro- 
tational symmetry involves a gapless Nambu-Goldstone mode 
P31l . while here we deal with a gapped system. 

We believe that this example constitutes a paradigm of 
mean field symmetry C/B in the course of adiabatic evolu- 
tion of a many-body system. The character of the strongly 
correlated states depend on the nature and character of the 
symmetry C/B, and the specific system - it is thus different 
in our case and in the case of Ueda et al. [45], or in the case 
of rotating lattice rings B51 . or in BEC in a tilted double well 
potential l47l . The last reference flTl reveals however, a fea- 
ture closely related to our results, which may have universal 
character: instability or chaotic behaviour of the system pre- 
dicted by the MF approach, is a signature of the existence of 
a strongly correlated state. 

This is quite an unexpected result, since at least for large 
systems, when TV goes to infinity, the MF theory in the regime 
of nucleation of the first vortex where the angular momen- 
tum of the systems changes from L = to N is believed 
to be correct. This belief has been in fact recently supported 
by rigorous results in Ref. [40 1, where it has been proved that 
the MF expressions for the total enerergy of the system coin- 
cides with the exact result for large N. Moreover, while the 
MF description of the GS is not correct for moderately big 
N at criticality, the single particle density is another example 
of a quantity correctly described by the MF approach. This 
last observation is supported by our result: we obtain that the 
density of the ground state at f2 c for increasing N becomes 
indistinguishable from that obtained for small N. This means 
that some macroscopic mean quantities, like the density or the 
total energy, are practically insensitive to the symmetry C/B 
process at criticality. 

The aim of the present work is to perform a deeper anal- 
ysis of the precursor state of the nucleation process, named 
\E' C from now on, using various techniques. First, we ana- 
lyze the dynamical process of increasing rotation frequency 
and show that, during the time evolution from an initial GS 
at a starting rotation frequency (where the angular momen- 
tum is L = 0) to the final one-vortex state where L = N, 
one must pass through a state that cannot be described by a 
MF order parameter. This state contains two macro-occupied 
modes of different parity of the single particle density ma- 
trix. Second, we characterize this non-condensed state and 
its properties, ask how does this state exhibit these proper- 
ties in measurement process. We note that similar questions 
have been posed in the context of appearance of the relative 
phase in the interference of two BECs in the course of mea- 
surements [48, 49, 50 1 . Inspired by the work of Javanainen et 
al. l48l we simulate the measurement process in a TOF exper- 
iment, for a one shot event, and for the acumulation of a large 
number of single shots. To be able to perform the measure- 
ment simulation of the \& c state assuming a large number of 
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particles, we use a two-mode model which provides a very ac- 
curate approximation of the exact state. Most of the times, the 
single shot events produce a single vortex located randomly 
along the x-axis whereas the accumulation of a large number 
of shots, reproduce the density. We discuss the interpretation 
of the outcomes and relate them with other strongly correlated 
states, discussed previously in various systems. 

Our paper is organized as follows: In Section II we present 
our model, and a brief repetition of previously obtained re- 
sults for rotating bosonic systems. In Section III we analyze 
the time evolution of the nucleation process. First we look at 
the possibility of an adiabatic evolution, and secondly, we dis- 
cuss a two-mode model for the GS at a critical frequency at 
which parity symmetry-breaking takes place at the MF level. 
In Section IV we study the energy spectrum as a function of 
the rotation frequency, in terms of the contributions of differ- 
ent L-subspaces in the GS and analyze the robustness of the 
^ c state. In Section V we describe measurement simulations 
and discuss their possible interpretations. Here the compari- 
son with a model cat-state (analogous to that predicted in Ref. 
|j46l ) is included. Finally in Section VI we present our con- 
clusions. 



and 

y 2 (z)=i^f>f (3 ) 

i 

where g — ^/8na/X z , a being the 3D scattering length, 
A 2 = \J H/Mlu z and Aj^ = ^/h/Mujj_. The dimensionless 
parameters g, A and B measure the strength of each term. 
We choose A^, huj± and uj± as units of length, energy and 
frequency respectively. It is worth mentioning that a simpler 
term that breaks the parity symmetry ~ B ^ i Xi is a center 
of mass excitation that would leave the internal structure un- 
changed revealing no new physics. In the second quantized 
formalism the Hamiltonian of the system projected onto the 
LLL in the rotating reference frame is described by: 



H = aL + PN + U + V 1 + V 2 = H + U+V 1 + V 2 , (4) 

where a — h(ui± — il) and j3 — fku±. HL , and N are 
the total z-component angular momentum and particle num- 
ber operators, respectively, Hq being the kinetic contribution. 
The contact interaction term is given by the operator 



II. MODEL AND BACKGROUND 

We assume a two-dimensional cloud of few condensed 
Bose atoms of mass M interacting via contact forces, confined 
in a symmetric parabolic trap with frequency ui±. Two extra 
perturbing potentials are also considered: One simulates a stir- 
ring laser that sets the system in rotation by a slight anisotropic 
deformation in the xy plane, rotating at angular frequency 
around the z-axis breaking the cilindrical symmetry. The sec- 
ond one, is a perturbation that breaks the parity symmetry, a 
symmetry that is otherwise preserved by the previous terms. 
The last one simulates possible second order contributions of 
the laser fields, and will help us in the analysis of the sys- 
tem. The rotation frequency 17 is strong enough to assume the 
lowest Landau level (LLL) regime [38], i.e., we consider that 
the kinetic energy within the LLL is given by H(oj± — il), the 
strength of the interaction, and both perturbations, are small 
compared with the separation between Landau levels given 
by h(ui± + f2) 1511 . It is implicit in our model that in the z 
direction, a strong parabolic trap of frequency uj z freezes the 
atomic motion producing an effective two-dimensional (2D) 
system. We model the contact interaction U, and the two per- 
turbing potentials V\ and V 2 by, 

N 

U(r t ,f}) = (h 2 g/M) ^Sifl-fj) ( 1 ) 

i<j 
N 

V 1 (r) = 2AMu J iJ2(x 2 l -y*) (2) 



U=^ Ul23i a Li a L 2 V"m 3 , (5) 

where the matrix elements read 

^1234 = (mi m 2 | U | m 3 m 4 ) 

g S mi+m2tm3+mi (mi + m 2 )! 
\ 2 ± ir VTOi!m 2 !TO3!m 4 ! 2"^+ m ^+ 1 

The operators a+ . and a mi create and annihilate a boson in 
a single particle eigenstate of H with angular momentum 
rrii, respectively. These eigenstates are taken as a basis to 
represent wavefunctions and operators in the second quan- 
tized formalism. We will refer to this set of functions as 
the Fock-Darwin (FD) functions restricted to the LLL, (with- 
out nodes in the radial direction) and will denote them as 
ip m = ^—^ (x + iy) m e~( x +y )/ 2A J- being m its angular 

momentum. The term V\ in Eq.4 is given by: 

Vi = — Y2 V m ( m - 1) a} n a m -2 

m 

+ v /(m + l)(m + 2) a) m a m+2 . (7) 
and the term V 2 by, 

V 2 = — ^2(^m(m - l)(m - 2) a^,a m _ 3 

m 

+ 3(m + aL a ™-i 

+ 3(m + 2)v / mTT a\ n a m+ i 

+ ^(m + l)(m + 2)(m + 3) al n a m+3 ) . (8) 
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In the absence of anisotropy (A = B = 0) the total angular 
momentum of the GS as a function of O shows sharp steps at 
critical values f2j i = 1, 2..., being f^i (Ox = <^_l — gN/&ir) 
the value at which the angular momentum of the system jumps 
from zero to L = N for all N. Above fij a plateau indicating 
constant angular momentum extends up to O2 where the sec- 
ond jump, not always to L = 2N takes place. From this value, 
a sequence of jumps and plateaux emerge up to the last pos- 
sible L-value given by L = N(N — 1), corresponding to the 
Laughlin state. The extension of the first plateau, from to 
O2 increases with N and at the same time, f2i decreases (this 
is a characteristic of the first plateau, in contrast, the next ones 
reduce drastically as N increases). From the expression of f2i 
it is evident that for a given value of g there is a maximum 
number of atoms compatible with our LLL assumption; for 
large N (N > 25 for g = 1) the functional relation between 
f2i and N cannot be linear any more, extra Landau levels must 
be taken into account. 

Eigenstates with more than one vortex can only exist lo- 
cally at Oj where several eigenfunctions with different angu- 
lar momentum, but degenerated in energy can be combined 
to generate vortex configurations in a spontaneous symmetry 
breaking mechanism. For instance, for N — 5 there is de- 
generacy between the states with L = 0, 2, 3, 4 and 5 at f2i, 
or L = 5 and 8 at ft 2 or L = 8, 10 and 12 at fi 3 . Namely, 
for small systems, the posibility to have vortex states is local- 
ized around discrete values of ft. However for larger N, on 
the one hand, the distance between the critical values of fij 
is drastically reduced [35 1, and on the other hand, the steps 
are softened due to the slight anisotropy that must be applied 
to accelerate the system; these two effects effectively provide 
the possibility to nucleate vortices in a continiuous way for all 

n>n x D21E5|. 

If anisotropy is included in the Hamiltonian, the steps are 
soften. For a fixed A (and B = 0) the softening is larger for 
small values of N in such a way that the steps may disappear. 
This effect can be reduced increasing g, namely, the increase 
of anisotropy reduces the effect of the interaction. To keep a 
well defined first step around the same it 1 for all N we con- 
sidered gN —constant in our calculations; in addition, this 
decrease of the interaction as N increases preserves the LLL 
condition [38 1. 

Having specified the Hamiltonian of the model, we per- 
form exact diagonalization. The isotropic Hamiltonian can 
be diagonalized in boxes of finite L-subspaces whereas, in the 
anisotropic case, several basis of different L-subspaces must 
be considered untill convergence is obtained, depending on 
the value of A and B. However, if B = since the term V\ 
in Eq. (7) can only mix L and L ± 2, the parity of the angular 
momentum is well defined and only subspaces with the same 
parity must be considered. 

To finish this section, we must introduce the single parti- 
cle density matrix (SPDM) operator and its eigenfunctions as 
tools to be used in the next sections. One of our aims is to 
analyse the states generated as ft grows in a dynamical pro- 
cess that follows the evolution of the GS from an initial value 



fli (smaller than a critical frequency fl c to be defined later) to 
a final value flf > ft c at the middle of the first plateau, where 
the first vortex has been nucleated. For these relatively low 
values of it and far from the critical frequency, the degree of 
condensation is high and the eigenfunction (a single particle 
wave function) of the SPDM that corresponds to the highest 
occupation plays the role of the order parameter (OP) of the 
condensate [52|. To obtain this OP we solve the eigenvalue 
equation for the SPDM, 

J d?n {l X?,?)r k {?) = n k ^ k {r), (9) 

where 

n^(f,P) = (GS I ¥(r)^(P)\GS), (10) 

with ^ being the field operator. If there exist a relevant eigen- 
value n\ >> n k for k — 2, 3, then 

v/mVitf) (li) 

plays the role of the OP of the system. The map of the lo- 
cal phase of this complex function gives precise information 
of the position of vortices [52|. The change of the phase by 
2irv around a point, where v is a positive interger signals 
the location of a vortex with v quanta of circulation. No- 
tice that m labels the single-particle angular momentum from 
rn = 0, 1, 2, ... of the FD functions, whereas k = 1, 2, ... is a 
label that distinguishes between the eigenstates of the SPDM. 
In Appendix A we show a detail derivation of the relationship 
between the operators that create FD functions: a m ^ and those 

that create eigenfunctions of the SPDM operator: bfj ■ 



III. DYNAMICS OF THE NUCLEATION OF THE FIRST 
VORTEX 

Adiabatic evolution 

First, we try to find out if the nucleation of the first vor- 
tex can be obtained as the final configuration of the adiabatic 
time evolution of the initial GS submitted to increasing il. To 
understand how the resulting state changes during the process, 
we perform several runs from ilo to ilf for increasing time in- 
tervals and compare the results with the sequence of station- 
ary solutions obtained from the diagonalization of the time- 
independent Hamiltonian, for the same range of O values. 

Assuming linear dependence in time of the rotation fre- 
quency as 

n(t)=ili + yt (12) 
where 7 is a constant, we solve the Schrodinger equation 

id t |$(t)) = H(t) |$(t)) (13) 
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FIG. 1: For TV = 6 and B — 0, time evolution of the mean value of 
the angular momentum, taking as the initial condition its GS-value 
at Qq = 0.4 (not shown in the figure). The value of At is the time 
used in the process to run linearly from f!o to fi/ = 0.85 in units 
of ljJ 1 . For N — 5, the inset shows < L > over Q; for 75 = 
once adiabaticity is fulfilled and close to adiabaticity for B — 0.003. 
At = 50000 in the last case. (A = 0.03 and g = 1 has been 
considered). 



for the initial exact solution \$>(t = 0)) at ilo- If we as- 
sume the expansion of the state in the Fock basis as | Q(t)) = 
Efc c k(t) \k) where \k) is the TV-particle Fock state with the 
well defined angular momentum, and project Eq.(13) on the 
state \j), we obtain the system of equations, 



id tCj (t) = Y,Ck(t)(j H(t) 



(14) 



which can be solved numerically. The bracket ( j H(t) k 
is expressed as 



H(t) 



L,(l-0(f))^ fe , (15) 



with a time-independent matrix which must be calculated only 
once, plus a time-dependent diagonal term. Hi = (3N + U + 
Vi + V'2\ Lj is here the angular momentum of the N-particle 
Fock state \j). 

We compare then the time evolution with the correspond- 
ing sequence of stationary states at instant rotation frequen- 
cies, and choose three different criteria of adiabaticity. The 
first one compares the profiles of the angular momentum as a 
function of ft. We consider that adiabaticity is obtained when 
the maximum difference between the curves is 0.03. In Fig. 1 
we plot for N = 6 (A = 0.03 and B = 0) the evolution of the 
angular momentum of the state that fulfills the condition of 
being the exact GS at fio = 0.4, up to ilf = 0.85 at the mid- 
dle of the first plateau. Different time intervals At have been 
considered untill adiabaticity is achieved for At = 30000; 
the black line corresponds to the stationary solutions. This 
means that from the stationary state at £Iq = 0.4 a slow evo- 
lution taking at least 28 seconds is necessary to nucleate the 



FIG. 2: Evolution of the weight (Pl) of the L-subspace components 
within $(t) before adiabaticity, for TV = 6 g = 1, A = 0.03 and 
5 = 0. 



first vortex at ilf evolving through stationary states (g = 1, 
A = 0.03 and u>± = 2ir x 175 Hz has been considered). The 
equilibrium state at Qo = 0.4 can be experimentally realized 
after relaxation, once the system is suddenly put in rotation at 
= 0.4 fl6l . The origin of the oscillations can be clarified 
by the analysis of the contributions of different L-subspaces in 
the expansion of |$(t)) as £1 evolves. Fig. 2 shows the weight 
of each subspace within it can be inferred that the inter- 
play between L = 6 and L = 4 produces the oscillations. We 
proved that this is a general result for all N, very fast evolu- 
tion is possible, keeping adiabaticity, from f^o U P to about f2i, 
whereas beyond this frequency, a sequence of quadrupolar ex- 
citations between L = N and N — 2 produce oscillations on 
< L > lowering the speed. 

A similar behaviour is found for our second criterion based 
on the evolution of the expected value of the energy. In con- 
trast to the previous case however, the coincidence with the 
stationary values is obtained much faster, (about At = 5000 
for the GS), producing the incorrect impression of adiabatic- 
ity, incorrect since other characteristics of the system, as is the 
case of the angular momentum, are still not reproduced. 

Finally, the third option is to measure the overlap between 
the exact GS and |$(i)) at Qf and to consider that one gets 
adiabaticity when the overlap is larger than 0.99. For N = 6 
adiabaticity is fulfilled in 21 sec, compatible with the result 
obtained from the first criterion. 

For odd values of N the time evolution driven by a parity- 
preserving Hamiltonian (with B = 0) cannot carry up < L > 
from to N, and a slight perturbation, which breaks parity 
symmetry is necessary. In other words, the sequence of sta- 
tionary solutions of the parity invariant Hamiltonian develop a 
first order transition at about fii where even L-subspace con- 
tributions are substituted by the odd ones within the composi- 
tion of the GS. In the inset of Fig. 1 we show, for comparison, 
for N = 5 the time evolution for the case when adiabaticity 
is achieved for B = 0, and nearly achieved for B — 0.003 
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FIG. 3: Dependence on TV of 7 C , the critical value of the varia- 
tion of with time, necessary to arrive to adiabaticity. B = and 
B = 0.001 has been considered for the upper and lower curves re- 
spectively, and A — 0.01 and gN = 6 in both cases. The horizontal 
line marks the asymptotic value for large number of atoms. 

(At = 50000 in the last case). In the first case, the expected 
value of the angular momentum is around 4, and no vortex is 
nucleated, whereas in the second case a one-vortex state with 
< L >= 5 is obtained. 

Figure 3 shows j c defined as f2/ = £lo + 7 C A£ as a func- 
tion of TV. For fixed values of the initial and final frequencies 
taken for all cases at ilo = 0.4 and ilf = 0.85, we increase 
At until adiabaticity is fulfilled (using the third criterion pre- 
viously mentioned), and from it, we obtain j c . In both cases, 
with and without the parity-breaking term in the Hamiltonian, 
7 C converges, meaning that even for large number of particles, 
the process is possible at a finite time interval. 



The two mode state 

To analyze the evolution of the GS in more detail (we fo- 
cus on the case B = 0) we look at the eigenfunctions of the 
SPDM (see Eq. 9), and their occupations n%. They provide 
an alternative representation of multiparticle states and oper- 
ators by the substitution of the FD functions ip m by ip^. The 
result obtained is that during the whole time evolution from 
Qq to Qf, only two of these single-particle eigenstates, ipi 
and ip2 (the most occupied and the next one), with occupa- 
tions ni and ri2 respectively, play a role in the GS. Moreover, 
through the whole evolution (ni+712) /N > 0.9, and as TV in- 
creases, this joint occupation is even larger. In Fig. 4 the adia- 
batic change in time of the first two occupations is shown. 
This result strongly suggests the substitution of the GS by a 
two-mode state in which only ipi and ip2 are considered. We 
define as "critical" (f2 c ) the frequency where the coincidence 
of the occupations n\ and 712 takes place. The decrease of 
the anisotropy A reduces the width of the critical regime, but 
the peaks in the population of ipi and ip2 always touch each 
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FIG. 4: Evolution of the mean value of the occupations n\/N and 
712 /TV in the GS, as a function of Q,, once adiabaticity is fulfilled for 
TV = 12. The peaks touch at Q c = 0.7759. The addition of both 
contributions is also shown. The inset shows the results for different 
number of particles. As TV increases, the condensation far from Q c 
increases, and at the critical frequency, the value of ni,2 is closer to 
0.5 meaning that the two mode-model is indeed better. The graphs 
are horizontaly shifted for clarity . A — 0.03 and gN — 6 has been 
considered. 
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FIG. 5: Weight versus Q of the FD functions ipo and if 2 (po and P2) 
in the expansion of the most occupied wave function ipi, below fl c 
as well as their joint contribution. TV = 12, A = 0.03 and gN — 6 
has been considered. 



other. In contrast, if B ^ 0, n\ is stricktly greater than at 
fi c , and although these modes are still the dominant ones, the 
peaks do not touch at any frequency. The two-mode model 
works well, but the system becomes slightly more condensed. 
However, within the region around il c , the system lies always 
beyond the regime of applicability of MF theory |44|. In the 
MF framework, the mechanism that triggers the nucleation at 
il c is related with parity symmetry breaking in the order pa- 
rameter, and this manifests as dynamical instability, i.e., any 
slight perturbation of the MF ground state drives the system 
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FIG. 6: Analysis of the GS for N — 12 at the critical frequency in 
terms of the square of the scalar products, P n =\ {n : ipi ; N — 
n : V2|*o) j 2 - Fig. 6a is for B = and 6b for B = 0.001. 
| (ME | GS) |= 0.92 for B = and | (PB \ GS) \= 0.840 for 
B = 0.001. A = 0.03 and gN = 6 in both cases. 



of 0.2 ■ 




FIG. 7: P n (defined in the caption for Fig. 6) versus n for different 
number of atoms for (B = 0.001). The case N = 16 and N = 18 
coincide and the result converges for larger N. For JV = 18, | (PB | 
GS) |= 0.000. A = 0.03 and gN = 6 has been considered. 



far off equilibrium (within the MF dynamics). 

Another relevant information, is that only the three first 
LLL single particle states (for m = 0,1,2) have a significant 
weight in the expansion of ipi or ip2- Below il c , ipi is a com- 
bination of (fo and ip 2 , at f2 c it changes its nature to a state that 
contains only the (pi and remains as that up to f2/. The sec- 
ond most occupied state, ip2> developes the opposite changes 
in such a way that they interchange their composition at the 
critical frequency. As a consequence, the density of the GS 
shows predominantly, two symmetric vortices, produced by 
the combination of ipo and cp^, that move from the edge to the 
center (up to a non-zero shortest distance), as approaches 
fi c and a single centered vortex, produced by ip\ from f2 c to 
ilf. Fig. 5 shows the evolution of the weights ( p and p 2 ) 
of the FD functions in the expansion of ipi up to the critical 
frequency where they disappear. Close to £Iq, ipi is essentially 
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FIG. 8: Overlaps, | {St \ GS) | as a function of tt, of the GS 
over different states (St): the maximum entangled (ME), the fully 
condensed (cond), the cat state (cat) and the fragmented state (frag) 
for TV = 12, A = 0.03, B = and gN = 6. 



ipo (and ipi is essentially the GS) producing a fully condensed 
state. Just before £! c it has both components with significant 
weights and at il c these weights are substituted by the weight 
of the unique component (p% (not shown in the figure). 

Had we considered the combination of ip and cp 2 for the 
whole range from f^o to £1 larger than ui± (with the addition 
of a quartic potential), in the spirit of a MF approach, then 
the curves po and p 2 in Fig. 5 would cross each other, and 
would end with final values of P2 = 1 and po = at some 
f2, meaning that a single double-quantized vortex would be 
produced at the centre of the condensate, in agreement with 
the results obtained by Saito and Ueda l53ll . The evolution 
would not fulfill the condition of adiabaticity and the last state 
would be an excited state. 

In order to perform simulations of TOF experiments in Sec- 
tion V assuming a large number of particles to have good 
statistics, we combine the results from exact diagonalization 
with the use a two-mode model in the following way: at il c , 
we analyze the overlap between the exact GS and A^-body 
two-mode Fock states of the type, 



|m,n 2 ,o,o...) 



(16) 



where n\ + n 2 = N, which will be abbreviated as \ni,n 2 ). 
Here we must clarify a point concerning the definition of ipi 
and ip2- At il c and B = 0, n\ = n 2 and the "most occupied" 
state can be both of them. As long as we concentrate in the 
analysis of the GS at il c , we choose for ipi the expansion in 
terms of ip and ip 2 when B = and for ip 2 the one repre- 
sented by ipi. For a slight B ^ 0, n\ > n 2 , and the previous 
choice is the right one. 

Fig. 6a shows | < GS\n, N — n > | 2 as a function of n 
for B = 0. The result means that in the GS there is a nearly 
uniform distribution among the different components, (better 
as increases) and in addition, that \GS) is very close to 
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the maximally entangled (ME) state constructed from even n 
values, since the overlap | (ME | GS) \— 0.92 is indeed 
large, here \ME) is defined as 

\ME) = [\N, 0) + \N - 2, 2) + ... + |0, N)]/y/N/2+l. 

(17) 

Note that we are using here the concept of entanglement for 
identical particles corresponding to the mode entanglement 
(for various ways of defining entanglement in systems of iden- 
tical particles see ll54l l55l . This analytic approximation to 
the GS allows us to simulate experimental TOF measurements 
with an arbitrarily large number of atoms, where the only in- 
gredient supplied by the exact analysis are the coefficients of 
the expansions of and ip2 on the FD functions. We use it 
in Section V. 

Fig. 6b shows the same results for B^O. The main differ- 
ence between the two cases consists on the expansion of ipi 
and V>2 in the FD states. In the parity broken case, in both ex- 
pansions, <po> f>\ an d f2 are significant. The decrease of the 
columns with n is exponential and can be accurately adjusted 
by an analytical function. We define as \PB) (PB for parity 
broken) the expansion on \n, N — n) states with the appropri- 
ate decreasing coefficients. It must be pointed out that the sin- 
gle particle odd occupations must be zero for B — and very 
close to it. The variation of the weights for different number of 
particles converges as N increases, as shown in Fig. 7 where 
the results from N = 16 and N = 18 coincide. In Section 
V, we use the converged coefficient to simulate the measure- 
ments. For a given N the distribution of weights shown in Fig. 
6, is robust against changes in the anisotropy strength A. 

For comparison, we also analyzed, for B = 0, the overlap 
of the GS with other celebrated states. One is the "cat state" 
represented by the combination |iV, 0) + |0,iV) |46|, which 
means that as a result of a single shot, the system can only 
appear as a full condensed state of each type with a probability 
of 50%, similar phenomenology has recently been analyzed 
in Ref. iBTll related with optical vortex cat states. Other is 
the so called "fragmented state" with only one component in 
its expansion, given by \N/2, N/2) with a 50% of occupation 
for each mode [50], and the last case considered is the full 
condensed state. In Fig. 8, for TV = 12, we show the overlaps 
as a function of il for each case, the system evolves from a full 
condensed state at f^o to a quite condensed one at flf passing 
through a strongly correlated state, very similar to \ME) at 
the critical frequency. 



IV. ENERGY SPECTRUM 
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FIG. 9: Lowest contributions on the energy spectrum as a function of 
Q, for N — 6 and gN — 6. In Fig. 8a we consider A — B = 0, in 8b 
A = 0.03 and B — and in 8c A = and B = 0.01 respectively. 
The arrows mark the value of fii. 




In Section III we addressed the question related to the pos- 
sibility of adiabatic evolution from ilo to ilf in a finite At, 
when non-zero values of A and B are included in the Hamilto- 
nian. The conclusion was that as long as B is small compared 
with A and g, even for a large number of particles, At remains 
finite. In this Section, in contrast, we analyze the energy gap 



between the GS and the first excited state as a necessary ingre- 
dient, aside from the critical 7 C already treated in Section III, 
to decide if the adiabatic criterium is fulfilled. We look for the 
possibility of excitations of different multipolarity. Firstly, we 
analyze the energy gap between the GS and the lowest excited 
states during the evolution in terms of the contribution of the 
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FIG. 10: Weight of the L-subspaces in the composition of the GS 
(fig. 9a) and the first excited state (Fig. 9b) versus Q. N = 10, 
gN = 6, A — 0.03 and B — 0.005 has been considered. 



L-subspaces in their composition, and see how this analysis 
depends on A, B and N. Next we concentrate on the robust- 
ness of the GS of the parity-invariant Hamiltonian at the crit- 
ical frequency, against occasional external perturbations, i.e., 
we inquire about the experimental feasibility of obtaining the 
$ c state in the laboratory. This question is suggested by possi- 
ble applications to quantum information, where it is necessary 
to manipulate and control the system. 

For the symmetric Hamiltonian (A = B = 0), as was men- 
tioned in Section II, at the critical frequency (f2 c = Oi = 




FIG. 11: Minimal gap (G m ) between the GS and the first excited 
state as a function of -/V. In the upper curve only even values of L 
has been considered (2x2), whereas for the lower curve, all L are 
included (1 x 1). A = 0.03, B — and gN = 6 are considered. 



uj± — gN/8n in this case), N stationary states are energy de- 
generated (L = and L = 2, .., AT) as expected from the 
analytical expresion of their energies given by, 



4tt 



L 
2 



1) + (lj±-Q)L + Nuj a 



= L{uj ± -SI-& + ^(N - 1) + Nu ± ,(18) 

which becomes L independent when il = f2]_. This degener- 
acy is lifted by the introduction of non-zero A or B or both. 
Fig. 9 shows for N — 10 the evolution of the energy spectrum 
for A = B = 0, A ^ and B = 0, and A = and B ^ 
respectively. In the second case, the energies are grouped in 
pairs of different parity, as was previously obtained by Parke 
et al. (ref), namely, for some values of il, the lowest excited 
state is nearly two-fold degenerated being however, well sep- 
arated from the GS, for the values of TV considered (N < 20). 
As O changes, some crossings and anti-crossings take place 
in such a way that the minimum gap between the GS and the 
first excitation (defined as G m ) may imply jumps from even 
to even (which is the case of N < 6), or even to odd (for 
N > 6). 

In Figure. 10 we show, for A^ = 10, the evolution of 
the contributions of different L-subspaces in the GS (Figure 
10a) and in the first excited state (Figure 10b). We consider 
A = 0.03, and a relatively large value of B — 0.005 to em- 
phasize the broken-parity effect. Since both symmetries are 
broken, all L-subspaces have non-zero contributions, however 
there is a remarkable difference between the regions below 
and above f^i = 0.76. For < Hi only even values of L 
are significant in spite of the fact that the Hamiltonian breaks 
parity symmetry, and the presence of external fluctuations that 
would produce monopolar (L ± 1) or octupolar (L ± 3) exci- 
tations is irrelevant, since these excitations are energetically 



blocked; the most probable process is a quadrupolar excita- 
tion from L = to L — 2. In contrast, for £1 > £!i all even 
and odd values of L play a role; in particular, far from £1%, the 
change from L = 10 to L = 9 in the most probable scenario, 
which means that the possible non-adiabaticity of the evolu- 
tion would be dominated by a braking process (i.e., lowering 
the speed), where in an effective way, one of the atoms jumps 
from the condensed to the thermal phase, ceasing its contri- 
bution to the total angular momentum of the system. Had we 
suppressed B in Fig. 10, we would obtain for the first excited 
state, separated regions with only even or only odd contribu- 
tions (due to crossings in the spectrum) and in particular at 
the critical frequency, the jump would be from even to even 
for N < 6 and even to odd for N > 6. 

Finally, in Fig. 1 1 we show the minimal gap G m as a func- 
tion of N obtained from the parity-invariant Hamiltonian. We 
distinguish between two cases: in the upper curve (2 x 2) only 
even values of L are considered, i.e., is the gap for quadrupo- 
lar excitations, in contrast, in the lower branch all i-values are 
considered (lxl) and represents the minimal gap that must be 
overcome by any perturbation of arbitrary multipolarity. If the 
system in protected against parity-breaking perturbations, the 
adiabaticity of the process is guaranteed (upper curve) since 
G m tends to a constant for large N and simultaneously 7 C de- 
creases, in such a way that the adiabatic criterium given by 
j c N/G m <C 1 is fulfilled. However, if parity breaking ex- 
citations can occur and the number of particles is large, the 
adiabatic evolution is practically impossible. 

It is worth noticing that we have been dealing with two 
different definitions of the critical frequency il c , one is the 
frequency where the two single particle occupations equal- 
ize (ni = 712), and the other one is the frequency where the 
minimum gap between the GS and the first excited state takes 
place; within our numerical precision both definitions are the 
same. 



V. SIMULATION OF MEASUREMENT 



Once we have an accurate representation of the GS at the 
critical frequency for both cases, with and without parity sym- 
metry breaking, we can simulate measurements during a TOF 
experiment for an arbitrary large number of atoms ll48l . We 
assume balistic expansion, a hypothesis compatible with our 
LLL condition ll56ll . We proceed as follows: after solving the 
eigenvalue equation for the SPDM of the exact GS (see Eqs. 
9 and 10), we get ipi and ipz that define the two modes. Using 
Eq. 17 we obtain two different realizations of the two-mode 
model, defined as \P0) for B = and \PB) (PB for parity- 
broken) for B / 0, both of them approximate the GS. Then 
the density distribution of this GS determines the position of 
the first atom using the following algorithm. A randomly gen- 
erated position r is accepted if p(r)/p max is larger than an- 
other randomly generated number u, < u < 1 or rejected 
otherwise; p ma x being the maximum of the density. This first 
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FIG. 14: A single shot on the \PB) state (PB for parity-broken) 
and its N-th correlation function for N = 1000, A = 0.03 and 
B = 0.001. 

step ends up with the detected position of the first atom, let's 
call it f\. Next we consider the pair correlation function (PCF) 
given by p (2) (fi,r) =< GS \ ¥{r^¥{f)^{f)^{r-s) \ 
GS > ( being 3?(r) = J^i <PlQl me ^sld operator ) and gener- 
ate the second position in the same way. By repetition of the 
procedure iV-times, finally we get the set of positions of all 
the N atoms. This sequential algorithm simulates a TOF mea- 
surement of a single shot. The procedure is the same in both 
cases, with or without parity broken symmetry however, the 
main difference is that in the first case (B ^ 0), the weights 
of the components \n, N — n) in \PB) become negilible for 
n > 18, which simplifies the numerical calculation. 

Fig. 12 shows a set of 4 shots for N = 10000 atoms with 
B = 0; here the spots represent the positions of the N atoms, 
whereas in Fig. 13 we show the iV-correlation function de- 
fined as, 

P {N \f 1 ,f 2 ,...,? N - 1 ,r) (19) 
related to the 4 shots of Fig. 12 respectively. The N- 
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FIG. 15: Shots on the Cat State. 



correlation gives the probability to find the last atom at r, 
once the other N — 1 have fixed positions. For such a large 
number of atoms we obtain, as a general result, that the j- 
correlation function is indistinguishable from the j + 1 one, 
if j > J ~ 100, namely the correlation is important for the 
first positions, but ceases to be modified by the addition of 
more fixed positions, due to screening effects, somehow this 
number J is related with the locality of the correlation. From 
the sample of 4 single shots shown in Fig. 12, it is clear that 
a vortex is produced at random places along the X-axis (this 
axis is preferred due to the particular form of the anisotropic 
term considered), with sometimes a slight manifestation of a 
second vortex, as is the case of the first shot. At odds with 
this result, in the case where parity symmetry has been bro- 
ken (i.e., the last degree of freedom), the GS is projected onto 
the most probable option and the vortex always appears at the 
same place (see Fig. 14), at a negative value of x in our case. 
This result is again due to our particular choose of the parity- 
broken term in the Hamiltonian, namely, the vortex would be 



A,r 2 N 



X 



(20) 



_L 



had been chosen. 

We interpret the results in the following way: each realiza- 
tion of the system in a single shot (in the case B = for ex- 
ample) is produced by two ingredients, one is the "intrinsic" 
nature of the GS which determines the density distributions 
of the j < N successive correlation functions, (or equiva- 
lently the density of systems with N — j + 1 particles), and 
the other is the particular measurement procedure, in our case 
by a one-by-one detection of the atoms. In other words, the 
measurement modifies the system and the pictures shown are a 
combination of the two factors. Different type of experiments 
detecting particle positions, would produce different pictures, 
coming however from the same GS. However, the differences 
disappear in the averaged picture given by the density, the ex- 
perimental mechanism does not modify the mean properties 
of the system. 

According to the result demonstrated in Appendix B, the 
accumulation of a large number of shots for a macroscopic 
system, reproduces the density. However, a non expected re- 
sult is that the density of a system of N = 10000 particles is 
indistinguishable from that of a reduced N (N < 20), mean- 
ing that some mean properties of macroscopic systems are 
well captured by the exact results from mesoscopic systems 
making unnecessary the extrapolation to the thermodinamic 
limit. This density contains two vortices symmetricaly posi- 
tioned along the x-axis. We want to stress that at the level 
of the exact GS we proved that the holes that appear in the 
density are real single quantized vortices as the phase of ifii 
changes by 2ir around them. If the analysis of the nucleation 
process would had only the density of the sequence of station- 
ary states for increasing fi from fin to fi c , as a unique source 
of information, the conclusion would be that the nucleation of 
the first vortex is preceded by the presence of two symmet- 
rically positioned vortices that move to the center. However, 
this is nothing else than a possibility within the multiple real- 
izations of the experimental performance. 

As a test of our procedure to simulate the TOF experiment, 
we considered a cat state, artificially created from the |P0) 
supressing all the contributions shown in Fig. 6a except those 
from \Q,N) and |iV, 0). The result is shown in Fig. 15, which 
corresponds to the densities of tp\ and ifa with one or two 
vortices as expected. We obtained only two possibilities, as 
the system is fully occupying the single particle ipi or ip2 and 
no partial occupations are possible. The positions of the non- 
centered vortices is given by x — (s/2 \ % I) 1 / 2 ~ 1.5 where 
tpi = atpo + /3(f2, in terms of the FD states. 

Some experience in the analysis of many -body systems sug- 
gests that complementary information can be obtained by the 
inspection of the PCF ||36l . However, the success of this ex- 
cercise may be strongly different for different type of states. 
Since the meaning is the probability of finding an atom at f 
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when another one is placed at fi, a circular symmetric re- 
sult means that the system is not correlated, whereas a hole 
at fi and ordered peaks in special positions, means that the 
position of the atoms is strongly correlated; this is, for exam- 
ple, the case of the Laughlin state, solution of the symmetric 
Hamiltonian with L = N(N — 1) in the region of strong rota- 
tion. This is a non degenerated GS, and as a consequence, its 
dendity is circular symmetric preserving the symmetry of the 
Hamiltonian. However, its PCF reveals a strongly correlated 
Wigner type structure of N peaks. The reason for the discrep- 
ancy between the density and the PCF is that the Laughlin 
state in a number state with the phase completely undefined 
(the opposite case of a condensate which is a phase state), it is 
a linear combination of all possible orientations of the Wigner 
structure. A measurement that fixes the position of one atom 
projects the system in a particular orientation, revealing the 
Wigner structure. In contrast, the GS at il c is the solution 
of a Hamiltonian that contains a rotating symmetry-breaking 
term (A ^ 0) and has a fixed orientation. The PCF does not 
provide then extra unknown information if f[ is chosen at the 
maximum of the density. 

Finally, we want to mention a speculation. In the case of a 
system submitted to both types of anisotropy (A ^ and B ^ 
0), and considering in general both signs for the parameter B, 
the density of the system would contain two vortices in fixed 
symmetric positions on the .T-axis. This state is a candidate to 
experience tunneling between the two single vortex states in 
agreement with the picture raised by Parke et al. [37| as the 
precursor mechanism for the nucleation. 

VI. CONCLUSIONS 

We have analyzed in this paper vortex nucleation in mezo- 
scopic 2D Bose superfluid in a rotating trap. The main 
ingredient of our work is that we have included a weakly 
anisotropic stirring potential, breaking thus explicitly the ax- 
ial (rotational) symmetry. The system we consider is well de- 
scribed by the mean field theory well below criticality (with an 
even condensate wave function), and well above the critical- 
ity, with the order parameter being the odd function. This phe- 
nomenon involves therefore a discrete parity symmetry break- 
ing. In the critical region the MF solutions exhibit dynamical 
instability. The main result of our paper is that the true many 
body state is a strongly correlated entangled state involving 
two macroscopically occupied modes (eigenvectors of the sin- 
gle particle density matrix). We have characterize this state in 
various aspects, which can be summarized in more details as 
follows: 

• The parity symmetry breaking at the critical frequency 
manifests itself as dynamical instability within the 
mean field framework. It does not prevent, however, 
the adiabaticity of the nucleation process. The increase 
of the parity conserving perturbation A notably reduces 
the necessary period of time to evolve from to 17/ 



(for A — 0.1, At ~ 1 s.). This conclusion remains valid 
even when the number of particles increases. However, 
a significant value of the parity breaking perturbation of 
B evidently acts against the adiabaticity. This perturba- 
tion leads to an exponential decrease of the energy gap 
from the GS to the first excited state. 

• The maximally entangled combination of ipi and ifj 2 of 
the two mode state, which is a fairly accurate represen- 
tation of the strongly correlated state at the critical fre- 
quency for B — 0, reveals a single vortex structure ran- 
domly located along the x-axis in a single shot measure- 
ment (with an additional small probability of a pair of, 
in general, non symmetrically located vortices). This is 
the result of the particular way of measurement mech- 
anism, that we consider here; a one-particle-followed- 
by-another-one detection. 

• The function p^ 2 ' (fi , r) with fi at one of the peaks of 
the density on the y-axis, does not reveal any hidden 
structure, due to the fact that the system has a fixed ori- 
entation, and the position of the vortex along the x-axis 
is smeared out by the integration over the positions of 
the other N — 2 atoms. This is an intrinsic property 
of this correlation function. In contrast, p N (r) breaks 
both symmetries, rotational and parity, producing the 
pictures shown in Fig. 13 typical of a projection mecha- 
nism implicit in a single measure |50|. 

• The state | ME) becomes a better representation of the 
exact GS as N increases. It is robust against changes in 
A for < A < 0.1. The state \PB) has zero contribu- 
tion of \n, N — n) for n > 18. 

• The mean properties of the system as those given by the 
total energy and possible by the density along the evo- 
lution in the whole range of variation of il considered, 
are insensitive to the symmetry broken mechanisme at 

n c . 

• Instability or chaotic behaviour of the system in a mean 
field calculation can be a signature of the existence of a 
strongly correlated state which description lies beyond 
the mean field framework. 



APPENDIX A 

Natural orbitals From the diagonalization of the SPDM, 
one obtains a new set of single particle wave functions (spwf), 
often refered as natural orbitals, and their corresponding ere- 
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ation and annihilation operators, 



3=0 
L 

hi = ^ftjfij 

3=0 



FD 



(21) 
(22) 



where are real numbers, hj is the operator that annihilates 
the FD with angular momentum j, j = 0, 1, L being L 
the largest single particle angular momentum in the GS \^o). 
We have sorted the spwf's in decreasing order of occupation 
(( rii )), in such a way that bi and b 2 create the most occupied 
single particle and the next one respectively. It is worth to 
notice that the subindex in b\ and ifii (i = 1,2, ...,L + 1) 
labels different SPDM eigenstates, whereas the subindex in 
ipf D , (i = 0, 1, L) means angular momentum. 

The representation of a state |W) in terms of the FD func- 
tions can be transformed into one in terms of the natural base 
{4>i} in the following way: 

Given a general state 



i*> = 5> fc |fc) 



FD 



(23) 



k=l 



where \k) FD is a iV-body state expressed in the FD base and 
dim is the dimension of the space considered, \k) FD can be 
expressed as, 



and m = 2, tpf PDM = c (p$ D + d ip% D , with d = -vT^c? 
and the second one is equal to the FD with angular momentum 
m = 1, if2 PDM = ff D - Then, each state can be expressed 
in the FD bases as 

|n, N-n,0, 0) SPDM - 1 2 1 1 , |0) = 

y/n\{N — n)\ 

E ,• ) A \0,N-n,0,-,0) FD - 




: d i \n-i,N -n,i,0, ...,0) 



FD 



(26) 



Projecting each of this states on the GS ob- 
tained from exact diagonalization \GS) = 
T,i , il ,...,i L a io,ii,-,ii.\ i o,ii,:;iL) FD , summing over 
all the states \n,N — n,0, --^spdM' n e even and finally 
multiplying by the normalization constant 1 / WN/2 + 1, we 
obtain the overlap expressed in the simple form 



[GS\ME)\ 



V N / 2 + 1 



N 



E E 

n=0 i=0 
nGeven 




C n 1 d l a n -i N- n i O,...,0 



(27) 



which measures the suitability of \ME) as an approximation 
oftheGS. 



N 



W)fd " 



I1«U|0) (24) 



where nit is the occupation of FD / in the state \k) FD , and 
Wkj is the angular momentum of each particle in the state 
\k) FD . From eq. 22 and eq. 24 we obtain for eq. 
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dim 

|*> = E 

*!=i 



N L+l 

, n e A,t 

lli=o 3=1 i=l 



(25) 



where we have used the ortogonality properties of the (3 
matrix: {fa}- 1 = {pij} T = {fa}. 



APPENDIX B 

The density The superposition of the data coming from a 
large number of single shots reproduces the density of a sys- 
tem in a state |\&), as can be shown in what follows. 

The probability to find a particle at the position r after k 
particles have been detected, is given by the function (k < 
N — 1, N the number of particles), 



,(*+!)/ 



(r 1 ,r 2 ,- - ,r k ,r) = (r^ (r 2 ) ■ ■ (r k )- 



(28) 



#t( r )£( r )#( rk )...#(r 2 )#( ri ) 



Overlap The ME state is a combination of states 
\n, N — n,0, ...) SPDM , n e even, with the same weight, 
where n is the occupation of the most occupied spwf and 
N—n the occupation of the next one. At £! c the first natural or- 
bit is a combination of spwf with angular momentuma m = 



where the expected value is in the state ^ and ^(rt) = 
Si^i^k)^ is the field operator. Using the commutation 
relations of the creation and annihilation operators and the 
ortogonormalization of the set {(pi(rk)} we can deduce the 
general relations 
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/ / ■■■/ 



/ p^{v 1 ,v)dv 1 - (N-l) p(r) 
/ / p( 3 )(ri,r 2 ,r)dnrfr 2 - (N-2)(N-l) p(r) 

p (JV) (ri,r 2 ,--- ,r N _i,r)dridr 2 ---dr N _i = (N - 1) ! p(r) . 



(29) 
(30) 

(31) 
(32) 



We could use any of them to recover the density, however the 
last is the one that fits our simulation. To model the experi- 
ment we have defined a grid in the xy-plane and we count the 
number of times that we detect a particle at each site of the 
grid. On the other hand, if we interprete the multiple integral 
of Eq. 12 as a multiple summation on a large number of dif- 
ferent configurations {ri,r 2 , ...,rN-i,r} on the discretized 
grid, keeping r fixed, then we can make the connection since 
the histogram obtained after a large number of shots, is noth- 
ing else than, aside of a constant number, the probability to 
find a particle at r when the other N — 1 have visited all the 
possible configurations, which is the meaning of the correla- 
tion function pW in the left hand side of Eq. 12. More than 
that, the summation over the p( N ^ functions is not arbitrary, 
it contains the information of the structure of the state >J/ as 
is the case in the simulation. It can be easily proved for ex- 
ample, in the case of the pair correlation function that can be 
rewritten as 



QUAGATUA, and Alexander von Humboldt Foundation Se- 
nior Research Prize. 



p< 2 >( ri ,r) = (Vb | *t (ri) ^t (r) ^ (r) ^ (ri) | 

= p(n) (Vi | *t( r )$( r ) | V \ (33) 



where 



\Vi) 



tt(ri) | Vb) 



(34) 



is a system with N — 1 particles and | Vb) is the initial state 
| VP). Or in other words, the probability to find a particle at r 
when other is located at ri depends on the probability to have 
a particle at ri in the state >!/. This complets our assertion. 
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